//
//  Simulation.cpp
//  phy_four
//
//  Created by Snir Gazit on 3/12/12.
//  Copyright (c) 2012 __MyCompanyName__. All rights reserved.
//

#include "Simulation.h"
#define NMODEL 2.

typedef boost::chrono::system_clock::time_point TimePoint;
typedef boost::chrono::duration<double, boost::ratio<3600>> Hour;
typedef boost::chrono::duration<double, boost::micro> MicroSec;

double calcW(double mu, double g, int m) {
	double c1, c2, c3, md = double(m), arg;
	arg = 0.25 * pow(mu, 2) / g;
	double f1 = 0.25 * (NMODEL + 2. * (1. + md));
	double f2 = 0.25 * (NMODEL + 2. * md);
	c1 = 0.25 * pow(g, -f1);
	c2 = -mu * gsl_sf_gamma(f1) * gsl_sf_hyperg_1F1(f1, 3. / 2., arg);
	c3 = sqrt(g) * gsl_sf_gamma(f2) * gsl_sf_hyperg_1F1(f2, 1. / 2., arg);
	return c1 * (c2 + c3);
}
void calcWeight(double muIn,  double gIn, int maxBond,
		ostream &outFile, BoltzWeightType &weight) {
	double c, n, mu = 2.* muIn , g = 4.* gIn;
	c = calcW(mu, g, 0);
	int m = 1;
	for (auto w = weight.begin(); w != weight.end(); w++) {
		n = calcW(mu, g, m++);
		*w = (n / c);
		c = n;
	}
}
// Worm Alg

double calcWRGauss(double mu, int m) {
	double fac = 1. + double(m);
	return pow(mu, -fac) * gsl_sf_gamma(fac);
}

void calcWeightGauss(double muIn, double g, int maxBond, ostream &outFile,
		BoltzWeightType &Rweight) {
	double mu = muIn;
	double n, c, w;
//	c = calcWRGauss(mu, 0);
	for (int i = 0; i < MAX_BOND; i++) {
//		n = calcWRGauss(mu, i);
		//	w = (n / c);
		//	c = n;
		Rweight[i] = (double(i + 1)) / mu;
	}
}
void WormAlg::init(QmcParameters &qmcParameters) {
	//Init Params
	qmcParameters_ = qmcParameters;
	string runNum = qmcParameters_.runNum.getVal();
	// OutFile
	string fileName = qmcParameters_.outPath.getVal() + "/out" + runNum
			+ ".txt";
	outFile_.open(fileName.c_str());
	outFile_ << "Starting N=2 Simulation, Run Number - " + runNum << endl;
	// Calc Weights
	//qmcParameters_.g.setVal(1. / (qmcParameters_.g.getVal()));

	// Init Lattice
	boost::array<int, DIM> sizes = qmcParameters_.Llist.getVal();
	lat_.init(sizes);
	lat_.setSectorType(ZSECTOR);
	// Init Rng
	eng_.seed(qmcParameters_.seed.getVal());
	// Init Moves
	//	createMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//	createMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	//  removeMoveXY_.init(&eng_, &lat_, qmcParameters_, XYWeight);
	//removeMoveIsing_.init(&eng_, &lat_, qmcParameters_, IsingWeight);
	createPairMove_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	shiftMove_.init(&eng_, &lat_, qmcParameters_, boltzWeight_);
	chooseLoopType_.init(&eng_, UniIntVariableParams(0, 1));
	setCoupeling(qmcParameters_.g.getVal());
	// Measurements
	measures_.init(&lat_, qmcParameters_, &boltzWeight_, &outFile_);
	// Read From File
	if (qmcParameters_.initFromFile.getVal())
		this->read();
	int totalNumberOfMeasure = (int) pow(
			(double) qmcParameters.binSize.getVal(),
			qmcParameters.numOfBins.getVal());
	//measures_.updateMinNumOfMeasure();
}
WormAlg::~WormAlg() {
	outFile_.close();
}

void WormAlg::thermalize() {
	wormLengthXY_ = acc();
	wormLength_ = 0;
	while ((ba::count(wormLengthXY_) < 10000))
		nextMove();
	double wormLength = ba::mean(wormLengthXY_);
	int closedRate = std::min(100,
			int(
					double(qmcParameters_.rateFactor.getVal())
							* double(lat_.getNumOfSites()) / wormLength));
	if (closedRate < 0)
		cout << "Error";
	measures_.updateClosedRate(closedRate);
	outFile_ << "Worm Length " << wormLength << " Closed Rate - " << closedRate
			<< endl;
}

bool WormAlg::run() {
	bool endTime, dumpFlag;
	namespace ba = boost::accumulators;
	ba::accumulator_set<double, ba::stats<ba::tag::mean> > wormLength;
	while (measures_.checkNotDone()) {
		nextMove();
		measures_.DoMeasurement(endTime, dumpFlag);
		if (dumpFlag) {
			this->dump(false);
		}
		if (endTime) {
			//this->dump(true);
			break;
		}
	}
	return !endTime;
}
void WormAlg::nextMove() {
	SectorType cSectorType = lat_.getSectorType();
	if (cSectorType == GXY) {
		shiftMove_.makeMove();
		wormLength_++;
		if (lat_.getSectorType() == ZSECTOR) {
			wormLengthXY_(wormLength_);
			wormLength_ = 0;
		}
	} else
		createPairMove_.makeMove();
}

void WormAlg::setCoupeling(double g) {
	if (qmcParameters_.simType.getVal() == "HEISENBERG")
		boltzWeight_.fill(1.);
	else if (qmcParameters_.simType.getVal() == "GAUSS") {
		calcWeightGauss(qmcParameters_.mu.getVal(), qmcParameters_.g.getVal(),
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzWeight_);
		qmcParameters_.U.setVal(0.);
	} else {
		calcWeight(qmcParameters_.mu.getVal(), g,
				qmcParameters_.maximalBonds.getVal(), outFile_, boltzWeight_);
	}
	shiftMove_.setCoupeling(boltzWeight_);
	createPairMove_.setCoupeling(boltzWeight_);
}

void WormAlg::simulatedAnnealing() {
	outFile_ << " Startting Simulated Anneling " << endl;
	double gHighT = qmcParameters_.gHighT.getVal();
	int numOfSMA = qmcParameters_.numOfGSA.getVal();
	double dg = (qmcParameters_.g.getVal() - gHighT) / double(numOfSMA - 1);
	double nLengthIsing, nLengthXY, nTotal;
	int numOfThermal = qmcParameters_.numOfThermal.getVal();
	vector<double> gSA;
	wormLength_ = 0;
	for (int i = 0; i < qmcParameters_.numOfGSA.getVal(); i++)
		gSA.push_back(gHighT + dg * double(i));
	for (auto it_g = gSA.begin(); it_g != gSA.end();) {
		double cg = (*it_g);
		setCoupeling(cg);
		outFile_ << "Current g -" << cg << endl;
		bool conv;
		int numIter = 0;
		wormLengthXY_ = acc();
		int i = 0;
		while (ba::count(wormLengthXY_) < qmcParameters_.SALoops.getVal())
			nextMove();
		//		if (lat_.sectorType_ != ZSECTOR)
		//			i++;

		int nLength = ba::mean(wormLengthXY_);
		nTotal = nLength;
		outFile_ << "Worm Length XY MA - " << nLength << " Count - "
					<< ba::count(wormLengthXY_) << endl;
		it_g++;
			//int fact=10000;

	}
	wormLengthXY_ = acc();
	acc lastTenAv;
	int count = 0,countten=0;
	for (int i = 0; i < qmcParameters_.SAAv.getVal(); i++) {
		while (ba::count(wormLengthXY_) < qmcParameters_.SALoops.getVal())
			nextMove();
		if (i>qmcParameters_.SAAv.getVal()/2)
			lastTenAv(ba::mean(wormLengthXY_));
		outFile_ << " Mean - " << ba::mean(wormLengthXY_) << endl;
		wormLengthXY_ = acc();
	}
	nLengthXY = ba::mean(lastTenAv);
	outFile_ << nLengthXY;
	nTotal = nLengthXY;
	double wormLength = nTotal;
	int closedRate = std::max(10,
			int(
					double(qmcParameters_.rateFactor.getVal())
							* double(lat_.getNumOfSites()) / wormLength));
	measures_.updateClosedRate(closedRate);
	outFile_ << endl << "Closed Rate - " << closedRate << endl;
}

void WormAlg::dump(bool final) {
	string outPath = qmcParameters_.outPath.getVal();
	measures_.dump(outPath, final);
	lat_.dump(outPath, final, &outFile_);
	eng_.dump(outPath);
	string fileName = outPath + "/sim_params.xml";
// Write to XML
	pugi::xml_document doc;
	pugi::xml_node rootNode = doc.append_child("SIM_PARAMS");
// Parameters
	pugi::xml_node node = rootNode.append_child("CLOSED_RATE");
	node.append_child(pugi::node_pcdata).set_value(
			(boost::lexical_cast<string>(measures_.getClosedRate())).c_str());
	doc.save_file(fileName.c_str());
//	doc.save(std::cout);
//	cout << endl;
}

void WormAlg::read() {
	outFile_ << "Starting Reading From File" << endl;
	string outPath = qmcParameters_.outPath.getVal();
	measures_.read(outPath);
	outFile_ << "Measures Loaded" << endl;
	lat_.read(outPath);
	int totalBonds;
	try {
		bool testConf = lat_.testConfiguration(totalBonds);
		outFile_ << "Total Bonds Calculated - " << totalBonds
				<< " total Bonds previous sim - " << lat_.getTotalBonds()
				<< endl;
		if (totalBonds != lat_.getTotalBonds())
			throw "Total Bonds after Read don't match";
		if (!testConf)
			throw "Bad Load of Configuration";
	}
	catch (string &e) {
		outFile_ << e;
		exit(1);
	}
	outFile_ << "Lattice Loaded" << endl;
	eng_.read(outPath);
	string fileName = outPath + "/sim_params.xml";
// Read From XML
	pugi::xml_document doc;
	try {
		pugi::xml_parse_result result = doc.load_file(fileName.c_str());
		if (!result)
			throw result;
	}
	catch (pugi::xml_parse_result &result) {
		std::cout << "Error description: " << result.description() << "\n";
		exit(1);
	}
	pugi::xml_node rootNode = doc.child("SIM_PARAMS");
	measures_.updateClosedRate(readXmlData<int>(rootNode, "CLOSED_RATE"));

}
